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We consider the thermal breakage of a tethered polymer chain of discrete segments coupled by 
Morse potentials under constant tensile stress. The chain dynamics at the onset of fracture is stud- 
ied analytically by Kramers-Langer multidimensional theory and by extensive Molecular Dynamics 
simulations in ID- and 3D-space. Comparison with simulation data in one- and three dimensions 
demonstrates that the Kramers-Langer theory provides good qualitative description of the process 
of bond-scission as caused by a collective unstable mode. We derive distributions of the probability 
for scission over the successive bonds along the chain which reveal the influence of chain ends on 
rupture in good agreement with theory. The breakage time distribution of an individual bond is 
found to follow an exponential law as predicted by theory. Special attention is focused on the re- 
combination (self-healing) of broken bonds. Theoretically derived expressions for the recombination 
time and distance distributions comply with MD observations and indicate that the energy barrier 
position crossing is not a good criterion for true rupture. It is shown that the fraction of self-healing 
bonds increases with rising temperature and friction. 

PACS numbers: 05.40.-a, 82.20.Uv, 82.20.Wt, 02.50.Ey, 



I. INTRODUCTION 



A great variety of problems both in material and basic science rely on the fundamental understanding of the 
intramolecular dynamics and kinetics of fragmentation (bond rupture) of linear macromolecular subject to a tensile 
force. Typical examples c omp rise material failure under stress [l], 0] , polymer rupture 0j7| . adhesion Q , friction 
mechanochemistry [T^, [Tlf , and biological applications of dynamical force microscopy [12i. [l3| . In particular, the 
problem of polymer rupture as a kinetic process has a longstanding history and dates back to the publications of 
Bueche 14 1 and Zhurkov et al. 15;]. In these papers the breaking of intramolecular bonds is treated as a thermally 
activated process and is described by Arrhenius's formula for the rate of bond scission. 

In recent years investigations have been complemented by ample application of computer experiments. Compre- 
hensive Molecular Dynamic (MD) simulations of ID chain fragmentation (at constant strain) have been carried out, 
whereby harmonic [Tg, [TtJ , Morse [l8l - [2"oT | , or Lennard- Jones [2ll - {24j chain models have been used. One of the principal 
questions to be answered is : "How long will it take for this system to break?" . A theoretical interpretation of MD- 
results [HI, H3] , based on an effectively one-particle model, has been suggested in terms of Kramers rate theory [25[ . 
Thus, the activation energy of the one-particle model Et, has been found to be close to the barrier height observed in 
the simulations, whereas the measured frequency of bond scission appears more than two orders of magnitudes smaller 
then the corresponding frequency, predicted by Kramers theory. The authors interpret this controversial result as a 
manifestation of the existence of collective modes which are missing in the one-particle Kramers theory. 

It has also been noticed that a truly irreversible break may occur, if bonds are stretched to lengths, considerably 
larger than the one corresponding to the barrier position. This indicates the possibility for bond recombination 
whereby the chain integrity is restored with some finite probability. 

The problem of polymer fragmentation has also been studied theoretically [2(| for the case of constant stress when 
a tethered chain of segments is subjected to a pulling force at the free chain end. The consideration has been based 
on a multidimensional version of the transition state theory (TST). Friction is then taken into account by coupling 
the polymer to a set of harmonic oscillators, simulating thus the presence of a thermostat. A comparison of the 
calculated breaking rate with the corresponding MD observation [21| shows again that the theoretical rate is about 
250 times larger. Also the role of the bond healing process has been discussed which helps to improve to some extent 
the agreement between theory and simulation. Nevertheless, despite the multidimensional nature of TST, it does not 
take into account properly the collective unstable mode development, which leads, in our opinion, to the essential 
overestimation of the breaking rate. 

The collectivity effect has been recently treated [26[ for constant strain and periodic boundary conditions (a ring 
polymer) on the basis of the multidimensional Kramers approach [23, Within this approach the development of 
a collective unstable mode and the effect of dissipation can be described consistently. It has been shown that in this 
case the effective break frequency is of the same order of magnitude as the one observed in the simulation. 

In the present work we develop this approach further for the case of a tethered Morse chain, consisting of N segments 
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and subjected to a constant tensile force / applied at its free end. We derive analytic expressions for the scission rate 
of the bonds, its distribution along the polymer chain, and its variation with changing temperature and dissipation. 
For comparison with computer experiment, we also perform extensive MD simulations in both one- ID, and three 
dimensions, 3D, and witness significant differences in the fragmentation behavior of the chain, despite the observed 
good agreement between theoretical predictions and simulation data. A major objective of the current study is the 
elucidation of the problem of bond recombination (self-healing) which has found little attention in literature so far. 
To this end we derive analytic expressions for the life times and extension distances of healing bonds and compare 
them to our MD results. 

The pap er is organized as follows. In Section QI] we give a sketch of the multidimensional Kramers-Langer escape 
theory [27, 28] . We also outline the problem of multiple points of exit from potential well [2!| [3(| which is necessary 
for treating bond rupture with respect to the consecutive number of each bond along the chain. In Section IIIII we 
present our model of one-dimensional Morse string of beads and consider the eigenvalue problem in the vicinity of the 
metastable minimum of the effective potential and at the barrier (saddle point) which is needed for the description 
of the unstable collective mode. Section llVl gives briefly details of our MD simulation and presents the main numeric 
results as well as their interpretation in the light of our theoretical approach. In Section [V] the healing process is 
discussed in terms of distributions of healing times and bond extensions. We give also the theoretical interpretation of 
this process based on the solution of the Kramers equation [31( using an inverted harmonic potential for representation 
of the barrier. We conclude in Section [VII and outline some future developments. 



II. KRAMERS-LANGER MULTIDIMENSIONAL ESCAPE THEORY 



A. Rate of escape 

The calculation of the rate of escape is based on the ./V-dimensional (in the total phase space Xi — qi, Xi+N = pi, 
i<N^) Fokker-Planck equation (3l| 



dP(x,t) 

at 



2N 
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(2.1) 



for the probability distribution function P(x, t) = P({xi},t). In eq. (I2.1[) the Hamiltonian has a general form 
YliLi Pi/2 m + The 2N x 2iV-matrix = — Aij , where is the matrix of Onsager coefficients and 

2N x 2N skew-symmetric matrix [32[ 
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(2.2) 



where 1 is the N x N unit matrix and is the N x N zero matrix. The eg. (12.11) can be seen 



as a continuity equation, dP/dt = — Y^i=i d/dxi Ji, where the probability current is given by Ji 



- E'fi M l3 {dH/dXjPfr, t) + T9P(x, t)/dXj). 

It is assumed [27[ that there is a metastable state {xf} which is separated with a barrier from another stable 
state. The coordinates at the saddle point which separates these two states is denoted by {xf}. The escape from the 
metastable minima is a comparatively rare event, so that one can treat the process close to {xf } as a stationary one, 
i.e. Y^j=i dJj/dxj = 0. Thus one gets 



2N 

E 



^2E%(x k -x s k )+T 
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dxj 



P(x) = 



(2.3) 



where also the harmonic approximation around {xf } has been used, i.e., the Hamiltonian reads 



2N 

H({ Xl }) = E S + £ E* k ( Xj - xf){x k - xf) 



j,k=i 



(2.4) 



where E s = H{{xf}) and the Hessian matrix Ef- = d 2 H / dxidxj at {x^ = {xf}. 
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One should impose the following boundary conditions. Near the metastable state {xf} the distribution function 
P(x) is the equilibrium one, i.e. 

P(x) = P cq (x) = Z A Y exp(-/3H) at x ~ {xf } (2.5) 

On the other hand, all states around the stable minimum (which is far beyond the top of the barrier!) are removed 
by a sink, i.e., 

P(x) ~ at { Xl } far beyond {xf } (2.6) 

With a transformation to new coordinates £ n = Y^=i Dni(xi — xf ), which are principal-axis coordinates for the 
Hamiltonian given by eq. (|2.4p . one obtains 



2N 

H(Z)=E S + -Y, + (2.7) 



71=1 



where the vector £ = {£„}. The matrix D n i is orthogonal one, i.e. D T = D^ 1 and {A„} are the eigenvalues of the 
matrix Since {xf } is a saddle point one of the A's, say Ai, is negative. The standard trick would be to look for 
the solution in the form P(£) = W(£)P cq (£) , where W(£) is a new function. Thus in the ^-coordinate system the 
steady - state Fokker Planck equation for W(£) takes the form 

|^ /- d 2 W ~ dW\ 

V r nfe — —-p\ n ^ n M nk — =0 (2.8) 



.fe=i 



where M nk = £\ j D ni M i:j D k j = T nk - A nk . 

In the same manner as for the one-dimensional Kramers problem [25[ one could claim that the function 
W(£) depends only on a linear combination of all £„, i.e. W({£ k }) — F(u), where u = J2' n U n £,n ■ The 
prime in this expression indicates that we omit all n's for which A„ = 0. The resulting equation reads: 

k (^nk U n U k ) d 2 F / du 2 — /3(A„ £ n M nk U k )dF/du = 0. As in the one-dimensional Kramers problem [25[ the 
coefficient in front of dF/du is a linear function of u, i.e. ^2 n k A„ (, n M nk U k ) — nu — k^2' n U n £ n . As a result 

A„ s ^M nk U k =KU n (2.9) 

k 

i.e., the coefficients U n are solutions of the eigenvalue problem eq. (|2.9p . Thus, 

, d 2 F(u) dF(u) 
T\ + —Y-u—^=Q (2.10) 
du du 



where 



A+ = - Y U n T nk U k (2.11) 



K 

n,k 



There is a simple physical interpretation of the eigenvalue problem given by eq. (|2.9p [27j . Namely, the equation of 
motion for the average value, (£ n (*)) = / E[?=i £nP({£j}, t) can be written as 



^<fn(<)>= JUdC j J({0,t) = ~^M nk ^JJdC—e-^ (2.12) 
With integration by parts and the harmonic approximation, eq. (|2.7[) . one arrives at the following expression 

§1 (£»(*)) = - E Mnk \k (&(*)) (2-13) 

The unstable solution of this equation (which describes the decay of a metastable state) is given by a negative 
eigenvalue k, namely 

(Ut))=X n e- Kt (2.14) 
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Substitution of eq. (|2.14[) in eq. (|2.13l) leads to J2 k M n k \k X^ — nX n or 

A„ ^nk A * X k = K\ n X n (2.15) 
k 

This equation is identical to eigenvalue problem, eq. (I2.9[) . provided that U n = X n X n . There is a negative eigenvalue k 
and as a result a negative A+ (see eq. ()2.11|) ) which corresponds to the unstable mode. This clear physical interpretation 
justifies the linear combination ansatz which has been used upon the derivation of eq. (|2.10p . We will show in Sec. 
IV devoted to MD-simulation results that the law given by eq. (|2.14l) actually holds for the breaking bonds. 
With the negative eigenvalue, A+ < 0, the solution of eq. (|2.10[) takes the form 



In eq. (|2.16[) we take into account the boundary conditions, eqs. (|2.5p and (|2.6I) which require F(u — > — oo) = 1 (at the 
metastable well) and F(u — > oo) — (around the stable minimum). 

Consider the rate of the metastable state decay. The main quantity which should be used for this purpose is the 
the probability current 

J « ^ -IT E ^n k U k e-P H (2.17) 

k 

To obtain the total probability flux over the barrier one should integrate the current, eq. (|2.17|) . over the hypersurface 
u = ^ Un^n — containing the saddle point. The resulting flux reads 

. 27V 2N 

J = / II d ^ S( -^ £ U n J n {{^}) (2.18) 

i=l n=l 

The calculation of this integral is given in Appendix of ref. 28] . The calculation yields 



V 2 H /^TX 1 / 2 q-PE 



(2.19) 



2tt \ |Ai | / lJ- 2 \ A„ J Za 
In order to obtain the rate constant one must divide the flux over the popularion ua in the metastable well 

j n ^ q (fc}) = n (^h ( 2 - 2 °) 



n A = 



where we have used the harmonic approximation near the metastable state {xf}, i.e. H = E A + Y«Zi K^l/ 2 - 
Combining eq. (|2.19[) and eq. (|2.20l) . one finds for the rate constant, k — J/ua, the following result [27J 



2tt 



det(E A /27rT) 1 1/2 



|det(E s /27rT)| 



e~ pEb (2.21) 



where the activation barrier Eh = E s — E A . We recall that A„ are the eigenvalues (without zero- modes) of the Hessian 
matrix E s at the saddle point x s whereas A^ are the the eigenvalues of the Hessian matrix ~E A at the metastable 



point x A . 



B. The eigenvalue problem 

Consider in more detail the eigenvalue problem given by eq. (|2.15[) . i.e.y\ M n k A& X k = nX n . In the initial n-space 
this equation reads 

J2M nk El r X r = KX n (2.22) 

k.r 
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where My = — Afj and the friction Tij as well as the matrix A^j are given as 



/ ••• - \ 



7717 
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The 2N x 2N Hessian matrix and the 2iV-dimensional column-vector are 



(2.23) 
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Xi = 
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qN 
Pi 

P2 

\Pn I 



(2.24) 



After substitution of eqs. (|2.23p and (I2.24[) into eq. (|2.22p . and exclusion of momentums {pi}, one arrives at the 
eigenvalue problem 



N 



3=1 % ' 

N x N - matrix 

where the notation = d 2 V/dqidqj\s is used. The corresponding characteristic equation [33[ reads 

JV 

det [mK 2 Sij — mjnSij + Vy] = [mn 2 — rn^K + Afc] = 

k=l 



(2.25) 



(2.26) 



where {A&} is a set of eigenvalues of the Hessian . As mentioned before, only one eigenfunction, say Ai, is negative. 
This eigenvalue determines then the equation for the transmission factor k, i.e. 



TO 



The negative solution reads 



Eq. (I2.26|) has been discussed first in ref. [34 



7 2 , |Ai| , 7 



7 <0 
to 4 



(2.27) 



(2.28) 



C. First-passage-time approach 



The method of Section III AI is a quite general approach to evaluate the rate of escape from the metastable state. 
The alternative to the flux-over-population method involves the concept of mean first-passage time. For an arbitrary 
stochastic multi-dimensional process x(t) the mean first-passage time (MFPT) T\(x) is defined as the average time 
elapsed until the process starting out at point x leaves a prescribed domain which includes the point x A of the 
reactant state. is sometimes called domain of attraction of the metastable state, and is shown in Fig. [1] 

For the FPT investigation it is convenient to use the Fokker-Planck equation (for the conditional probability 
P(y, t; x, 0) of visiting the point y G fi at time t, provided it starts from x e il at t — 0) which includes the so called 
adjoint operator and has the form [3lj 



d_ 

Of 



P(y,t;x ) 0)=£+(x)P(y,t;x,0) 



(2.29) 
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FIG. 1: Exit from the phase space domain Q. Vectors denote the restoring force directions. The trajectory starts close to the 
metastable minimum and exits through a particular exit (saddle) point. 



where the adjoint operator 

Lt(xH£^)A +£ £A,(x)^ (2.30) 

In the present case the drift and diffusion coefficients are respectively: Ki(y) — — J2j MijdH/dyj and Dij — . In 
eq. (|2.30p e is responsible for the noise intensity (in our case e plays the role of the temperature) , so that at e — » 
the noise is weak. 

The n-th order moments of the first-passage time, T„(x), may be iteratively expressed in terms of the adjoint 
operator, [3l[ 

L'(x)T n (x) = -JiT n „i(x) , at x G Q 

T„(x) =0 , at xeffl (2.31) 

where the second equation simply means that the first passage time is zero, provided the trajectory starts at the 
separatrix dfl. The hierarchy given by eq. (|2.31[) should be supplemented by the initial condition To( x ) = 1 which is 
evident from the normalization condition for the FPT distribution. Therefore, the equation for the mean first passage 
time (MFPT), Ti(x), reads 

L f (x) Ti(x) = -1 , at x G O 

Ti(x) = , at x G dfl (2.32) 

On the other hand, for small noise, e — > 0, a trajectory starting within fl will typically first approach the attractor 
and stay within its neighborhood for a long time until an occasional fluctuation drives it to the separatrix dfl. 
Hence, MFPT Ti(x) assumes the same large value r everywhere in Q except for a thin layer along the boundary dfl. 
Under these condition the solution of the hierarchical equation, eq. (|2.31[) . can be obtained in the following for m [35j : 
Ti = t, T2 = 2t 2 , T3 = 2 ■ 3r 3 , . . . , T n = n! r". In result the probability density of the first passage times reads [35J 

p(t) = r^expi-t/r) (2.33) 

(we recall that r _1 t n exp(—t/r)dt = nW n ). The validity of this PDF has been shown in our MD-simulation (see 
Sec. IV). Finally, it can be proven [28|,[35[ that the Kramers escape rate k = (2r) _1 . 



D. Distribution of exit points 



Assume that there are a number of saddle points, b r where r — 1,2,.. . M, (referred to as exit pomts)which lie 
on the separatrix. We study then the exit events (see Fig. Q] where only two exit points are shown) . One may ask 
: "What is the probability to exit through a particular exit point h s G <9f2 irrespective of the time it takes?" This 
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problem has been discussed first by Matkowsky and Schuss [29j and later by Gardiner [30]. Here we give an explicit 
solution which may be compared to MD-simulations. 

The probability for escape through an exit point h s if the random trajectory starts at x € is defined as 

oo 

7r(b s ;x)= f dt'J2M*> S )Mh s ,t'\x,0) (2.34) 
o 

where Vi(h s ) is a component of the vector normal to the separatrix unit vector at h s pointing out of region fi. 
The flux Ji(h s , t'\x, 0) counts only trajectories which start at t = in x e fi and approach the exit point h s at 
the boundary (separatrix) at time moment t'. The integral over time in eq. (I2.34[) implies that this probability is 
calculated irrespective of the time needed for escape. 

The probability ^(b^jx) is governed by the backward stationary FPE, i.e. (see Sec. 5.4.2 in [30]) 

Lt( x ) 7r (b s ;x) = (2.35) 

The boundary conditions are: 7r(b s ; b s ) = 1 and n(h s ; x) = for any x E <9f2 if x ^ h s , i.e., 

vr(b s ;x) = <5 s (b s -x) for x e dfl (2.36) 

Using the same arguments as in Sec. Ill CI one can show that in the small noise limit, e — > 0, the probability 7r(b , x) 
remains the same everywhere inside f2 apart from a thin layer along the boundary dfl. Now we specify the general 
cq. (|2.35[) for the special case when the drift term has the form of potential (i.e., it is derivative of the Hamiltonian) . 
Then eq. (|2~35]l reads 



E 



, BHd „ d 2 



dxj dxi dxidxj 



7r(b 5 ;x)=0 (2.37) 



Close to the saddle point the Hamiltonian can be treated in the harmonic approximation, eq. (|2.4[) . and the drift 
velocity in eq. (|2.37|) becomes 

f) J-f 

K l = -J2 Mi~ = -J2 M t] E%{x k - bf) (2.38) 



Since the solution changes only within a thin layer along the boundary (in direction normal to the boundary 
layer!), one may introduce new local coordinates {z,y p } where z measures the distance from h s and {y p } is the set 
of tangential variables measuring the orientation around h s . The coordinates z — z(x) and y p — y P (x) are chosen so 
that 

Vz(u) = i/(u) 
i/(u)-Vy p (u) = 

z(b s ) = (2.39) 

where u g dft. The first equation in (j2 .39[) means that the coordinate z changes in direction of the t'-vector. The 
second equation implies that the coordinates {y p } are parallel to 90. In terms of new variables 

V ^=^^+E v ^( x )^ (2 - 40) 

and 

+ ViVj-zfx) 97 + E V « V ^p( x ) ( 2 - 41 ) 

One should keep in mind that for the exit event at e — > only one coordinate is relevant, namely, the one traversing 
the saddle point in direction of v vector. This is our new z-coordinate, which parameterizes the displacement from 
the saddle point as follows 

z = ^fep (2.42) 
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Substituting eqs. (|2.40[) - (12.421) in eqs. (|2.37|) and (|2.38l) . and keeping only the lowest order in s, leads to 



where 



<9 2 7r dn 



, 1 r ™7 



(2.43) 



(2.44) 



In (!2.44[) one has used J2i j v i Mij Vj = J2i ? an d taken into account J7 j v i ^ij v j = TO 7 (Appendix [5} as 



well as 



1,3 



k= s }^ j i>j MijEf k v k 



(2.45) 



The solution of eq. (|2.43[) has the form 



ir(b s ; u, p) = 5(u - h s ) + [Coo - S(u - h s )} 



7r|A + | 



1/2 



exp 



2 A-i 



(2.46) 



where we took into account the boundary condition 7r(b s ;u, 0) = <5(u — h s ) (vector u e 90) and at p —> — oo 
7r(b s ) = Coo, i-e., well inside the region f2 the solution is a constant as this should be for a weak noise. In order to 
fix the constant Coo one may multiply eq. (|2.37[) by p st (x) and after integrating over Q (using also the integration by 
parts and the Gauss theorem) derive the surface integral 



j dS I - Pst J2 ViMijVjH n + e 



= 



Given that p st = Z A exp[— H/e], the 1-st and 3-rd terms in eq. (I2.47P cancel each other and therefore 

J dS p st ViMijVjn = 
an l ' J 

The gradient Vj7r, calculated from ea. (|2.46l) at the boundary (i.e. at p = 0), reads 

1/2 



(2.47) 



(2.48) 



VjTii^o = Vj [Coo - S(u - b b )} 



7T A, 



Substitution of eq. (|2.49[) into eq. (|2.48l) leads finally to the result for the exit probability 



Tr(b b ) = C 



-H(b s )/e 



M 



(2.49) 



(2.50) 



£ v/MV)e- ff ( br )/ e 



r=l 



In the last equality in eq. (|2.50[) the surface integral is replaced by a sum over all saddle (exit) points. This is possible 
because all eigenvalues of the Hessian ViV :( -iJ(u)| u= b>' are positive within the boundary hypersurface dQ and hence 
H(u) has minimums at u = b r (where r — 1,2, ... , M). The exit probability given by eq. (12.50)) will be compared in 
Section ITVl with our MD-simulation results. 

Finally, we stress that the choice of the coordinate system given by eq. (|2 .39[) so that only z-direction is physically 
relevant is in complete agreement with the Kramers-Langer approach where the linear combination (see the paragraph 
after eq. (12.81) ) plays the role of the relevant coordinate. Indeed, the first relationship from eq. (|2.39[) can be formally 
solved as 



*( x ) = v A x i ~ b 



(2.51) 



which after transformation to the principal-axis coordinates, i.e., £j = Y2j ^ij( x j ~ bf) leads to z({£j}) = y*, Uj£j 
(where Uj = Djk^k)- Moreover, on the separatrix y^. Uj£j = and one recovers the Kramers-Langer choice of a 
relevant coordinate (see the paragraph after eq. (|2.8[) ~) as linear combination of all £j. 
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III. BREAKAGE OF AN ONE-DIMENSIONAL STRING OF BEADS 

Here we describe chain breakage by means of the Kramers approach, sketched in Section [TT1 

A. Model 



We consider a tethered one-dimensional string of N beads which experiences a tensile force / at the free end as 
depicted in Fig[21 Successive beads are joined by bonds, governed by the Morse potential, Um(u) — D(l — e~ ay ) 2 , 
where D and a are parameters, measuring the bond strength and elasticity. The total potential energy is 



N 



V{{xi\) = ^ Um{x % - Xi-i) - fx N 



(3.1) 



where we set xo — (see Fig. [2]). Upon change of variables, yi = Xi — Xi-i, one gets 



JV 



JV 



V({Xi}) = pMiVn) ~ fVn] = Y, 



(3.2) 



so the combined one-bond potential then reads U(y) = D(l — c ay ) 2 — fy. 



^x Xj x 2 



-#^\^# 

x n x n+ i 



l N-l 



FIG. 2: Schematic representation of a tethered string of beads, subject to a pulling force /. The corresponding coordinates 
are marked as xi, X2, ■ ■ ■ xn- A single "endangered" bond (red) with length close to the distance of the one-bond potential 
maximum is located between the n-th and (n + l)-th beads and its spring constant K s may become negative. 



Direct analysis of the one-bond potential U (y) indicates that the positions of the (metastable) minimum y_ and of 
the maximum y + are given by 



y-,+ = - In 

a Ll±^l-/J 

where the dimensionless force / = 2/ / aD. The activation energy (barrier height) is given by 



(3.3) 



Eb = U{y+) — U{yJ) = D ■ 



i- Vi-/ 



i 



i-/ 



(3.4) 



One can easily verify that E\, decreases with /. Since the Kramers' theory implies 3> fc^T, the force / should not 
be too large. The characteristic frequencies at the minimum and maximum of U(y), ilf = (l/m)(d 2 U(y)/dy 2 )y = y_ 



and fl 2 



{l/m){d 2 U{y)/dy 2 ) y= 



v+ i 



are therefore 



n 2 - 

"1,2 — 



a 2 D 



1-/±(1-/) 



(3.5) 



Fig. [3^ illustrates the Morse potential as well as the modified one-bond Morse potential U(y) = D(l — e~ ay ) 2 — fy. 
In order to test the role of anharmonicity (recall that the Kramers-Langer theory uses only harmonic approximation) 
we have also tested in our MD-simulation the Double-Harmonic potential constructed piecewise as V(x) = (x— l) 2 — fx 
for x < 2 and V(x) — 2— (x — 3) 2 — fx for x > 2 and shown in Fig. [3t>. The corresponding barrier heights dependences 
on the pulling force are shown in inserts. 
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FIG. 3: Combined bond potential: (a) The pulling force gives rise to a metastable minimum and a barrier Ef,; the barrier hight 
Eb declines with growing force / as shown in the inset, (b) The effective Double-Harmonic potential constructed piecewise as 
V(x) = (x — l) 2 — fx for x < 2 and V(x) — 2 — (x — 3) 2 — fx for x > 2 for different values of the force /; the barrier hight Et 
declines with growing force / as shown in the inset. 



B. Eigenvalues close to the metastable minimum 



One can verify that the determinants of the Hessian matrix, entering eq. (|2.21[) . may be calculated exactly for our 
one-dimensional model and even the total eigenvalue problem may be solved analytically. 
The Hessian at the metastable minimum x A has a 2N x 2N block-matrix form 



/ 



E- 



d 2 H 



dx;dxj 






1/m 



\0 ■■■ 



• \ 

• 

• 

• 1/m J 



(3.6) 



where is the N x ./V-matrix of the potential energy second derivatives. It has the following tridiagonal structure 



/ 2 

-1 



V A 



1 

2 -1 
0-1 2-1 



-1 



V 



-1 
1/ 



(3.7) 



The correct diagonal and off-diagonal elements follow from the double derivation of V {{xi}) — U (xi) + U {x2 — x\) + 
■ ■ - + U(xn-i~ xn-2) + U(xn — xjv-i). For example, the upper-left element V$ = d 2 U{x\)/dx\-\-d 2 U(x2 — x 1 )/dx\ = 
2mVL\ whereas the bottom-right element reads V^ N = d 2 U{xM — XN~i)/dx% = m£l 2 . The only nonzero off-diagonal 



elements are 



V A 

v ii+l 



d U(x,i + i — Xi) / dxi+idxi = —mVL\. 



The block-matrix structure E^t- = ( q d ) ma kes it possible to calculate the determinant as [35 



det(E^ 



det A det D 



Thus, for the block-matrix in eq. 



one gets 



(3.8) 



det(E^) =[-) det(T^) 



JV 



(3-9) 



The calculation of the tridiagonal matrix is given in the Appendix [Bj Using eq. (|B5|) . one derives 



det(V A ) = {mVl 2 ) N and det(E 



2N 



(3.10) 



11 



The eigenvalue problem for the Hessian eq. (|3 . T[) reads V k A n 



Uk+l + U k -x 



Uk = 



(3.11) 



which should be supplemented by two boundary conditions, namely, uq — (tethered left end of the chain), and 
un = un+i (free chain end right). With the identity sin(fe + l)ipj + sin(fc — l)tpj — 2coscpj sm(ktpj) = 0, and 
comparing this with eq. (l3.11[) , one gets for the eigenvalues Xj = 2mf2f(l — cosfj). For the eigenfunctions one has 
Uk = sm(ktpj) with ifj being the mode factor which can be fixed by the free-end boundary condition tt/v = iijv+1- 
Eventually, one obtains ipj — (2j — 1)/(2N + l)w so that the eigenvalues read 



1 — cos 



[ 2j - 1 
V2iV + 1' 



(3.12) 



The corresponding eigenfunctions are 



CO 

= sin 



(2j - l)k 



2N + 1 



(3.13) 



C. The determinant of the Hessian matrix at the saddle point 

The chain breaks when at least one bond length comes close to the barrier position y+ , eq. (|3.3|) , of the one-bond 
potential U(y). The spring constant K s of this "endangered" bond is negative, i.e. 



a 2 D 



!-/+(!-/) 



(3.14) 



Such a bond, located between the n-th and (n + l)-th beads, is illustrated in Fig. [3] 
As before, (cf. cq. (|3.9p ). one has 



det(Ef ? ) 



N 



det(Vg) 



(3.15) 



where now the Hessian of the potential energy at the saddle point, V$ , has a more complicated structure than V^. 
Indeed, from the potential energy V = U (x 1 ) + U(x 2 — xi) + - ■ ■ + U(x n — x n -i) + U(x n+ i — x n ) + U(x n+2 — x n+ i) + ■ • • + 
U(xn-i — xm-2) + U{xm — CEjv-i) one recovers the following structure: Diagonal terms: V^f = 2mil\ fori ^ n,n+l,N, 
V-n+i.n+i — m (^i — ^|)> an d Vnn — m ^i] Off-diagonal non-zero terms: Vf i+1 = VA xi = —mill f° r * 7^ 71 > 



V ti = 

nn 

and V^ n+1 = V^ +l n = mil^. In result the tridiagonal Hessian matrix V/j reads 



( 2 1 

-1 2 -1 
0-12-1 



-1 1 




\ 



a a 
: 1 - , 



V 



-1 2 -1 
-1 1 / 



(3.16) 



(3.17) 



where a = il\/il\ < 1- Let us calculate first det(V^). The matrix can be considered as a block-matrix 



(3.18) 



where as arguments in (N, n) we keep the total matrix dimension, N, and the position of the "endangered" bond 
n. The n x n-block A, n x (jV — n)-block B : (N — n) x n-block C and (N — n) x (N — n)-block D in eq. (|3.18p are 
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given by 



.4 



( 2 1 

-1 2 -1 

0-1 2-1 



V 



0-12 -1 
-1 1 - a ) 



B 



••• 
a ■■■ 



(3.19) 



nx(N-n) 



c = 



••• a 
••• 



/ 1 - Q -1 

-1 2 -1 
0-1 2-1 



D = 



(N-n)xn 



\ 



0-1 2-1 

V o -11/ 



(3.20) 



(AT-7l)X(./V-7l) 



The block-matrix's determinant is given by 



det(V;f (N, n)) = (tnfl?) ff det(A) det( £> - CA^g ) 
The (N — n) x (N — n)-matrix F = D — CA^ 1 B can be readily calculated to 



(3.21) 



F(N - n) 



( Xn -1 

-1 2 -1 
0-1 2-1 



0-1 2-1 

V -11/ 



(3.22) 



where x„ = 1 — a[l + a(A 1 ) nn ], and (A 1 ) nn is the bottom-right element of the matrix A 1 . It is easy to verify 
that (A -1 ) ?m = n/(l — an), so one gets 



Xn = 1 - a 



an 



1 — an 



= 1 - 



1 — an 



(3.23) 



To calculate det(-F), the determinant is expanded in minors regarding the first row (see a similar expansion in Appendix 
dBj).This yields 

det(F(N - n)) = ( Xn - 1) . (3.24) 
On the other hand, by making use eq. (|B5|) . we have 

detL4(n)] = (1 - an) (3.25) 
Taking into account eqs. (|3.23[) . p.24j) and (|3.25[) in eq. (|3.21|1 . one obtains eventually 

det[V$ (N,n)] = -a{mn\) N . (3.26) 
The result given by eq. (|3.26p is only valid for 2 < n < N — 2. The cases for n — 0,1, N — 1 should be considered 
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separately. The Hessians in these cases look as follows 



V*(N,0) = mtt\ 



/ 1-Q -1 

-12-1 
0-1 2 



Vg(N,N-l) 



mill 



V 

/ 2 1 

-1 2 

-1 



V 



-1 



2 -1 
-1 1/ 



/ 1 — a a 

a 1 — a —1 



t$(JV,l) = mJ}? 







-1 



-1 

2 -1 



-1 



2 -1 
T 1 — a a 

a —a / 



0-1 2-1 
-1 1/ 



(3.27) 



Direct calculation which uses the Laplace's formula for determinants (see Appendix |Bj leads to the result 

det[V|(iV,0)] = det[^f(iV,l)] = det^f (N,N- 1)] = -a(mti{) N (3.28) 

Thus, comparing eq. (|3.28[) with eq. (|3.26p . one can see that the value of det[V^?(7V, n)] does not depend on the 
endangered bond index n. Taking into account eqs. (|3.15[) . (|3.26[) and (|3.28[) leads to the final result for the determinant 
of the Hessian matrix 



det[E s (iV,n)] = -aftj . 



(3.29) 



The determinant in eci. p.29[) is negative as it should. 

The ratio of the fluctuating determinants R(n) = [det E" 4 (A r )/| detE s (iV, n))] 1 / 2 which is involved in the general 
expression for the rate constant, eq. (|2.21|) . is given by 



R(n) 



(3.30) 



We emphasize that the ratio of the fluctuating determinants R(n) does not depends on n. The n-dependence of the 
total rate k is present in the k- factor (see eq. (|2.2ip ) which will be discussed in the next Section. 

D. The unstable mode 

In an unstable equilibrium configuration when an "endangered" bond has a negative spring constant there exist 
N — I stable modes and one unstable mode. One may find the eigenvalue A < and the eigenfunction Uk for the 
unstable mode. 

The eigenvalue problem for the Hessian, given by eq. (|3.17[) , reads Vj? r u r = Xuk- In detail this yields 



Ufc-i + Uk+i - I 2 - 



A 



to£1 2 



w«-i - au n+ i 



1-a 



for k 7^ n, n + 1 
for k — n 



-au n + u n+2 -1-a 



mil 2 



U n +1 



= for k = 



1 



(3.31) 
(3.32) 
(3.33) 



which should be supplemented by the boundary conditions: uq = (tethered end of the polymer) and un — un+i 
(free end of the chain). We recall that the index n denotes an "endangered" bond location between the n-th and 
(n + l)-th beads. 

Let us first find the solution for k < n. To this end we use the identity sinh[(fc + l)(p] + sinh[(fc — 1)99] — 
2 cosh(y>) smh(kip) = 0. Comparison of this identity with eq. (|3.31l) suggests that the eigenfunction which char- 
acterizes the distance from the unstable equilibrium position is given by 



Uk = — smh(k(p) 



(3.34) 
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with an eigenvalue 

A = 2mnf(l — coshtp) (3.35) 

where ip is the mode factor which will be fixed below. The solution eq. (|3.34p also meets the boundary condition 
u = 0. 

In order to find the solution at k > n + 1 we consider two identities 

sinh[(iV + 1 - k - l)<p] + sinh[(iV + 1 - k + l)cp] - 2 cosh(^) sinh[(A + 1 - k)y] = (3.36) 
cosh[(A + l-k - l)ip] + cosh[(A + l- k+l)tp]-2 cosh(y>) cosh[(A + 1 - k)ip] = (3.37) 

Again, comparison with eq. f|3.31[) suggests that there are two linearly independent, i.e., fundamental solutions 
( see, e.g., [Hj]), = sinh[(iV + 1 — k)ip] and = cosh[(iV + 1 — k)tp]. Thus the general solution reads: 
iik = Asmh[(N + 1 — k)ip] + B cosh[(N + 1 — k)ip] where A and B are some constants. The boundary condition 
un+i = ujv helps to express B in terms of A which gives 

u k = A |sinh[(iV + 1 - k)ip\ - coth (|) cosh[(iV + 1 - k)ip\ } (3.38) 

We extend the solutions, given by eq. p. 341) and eq. (|3.38[) . up to k — n and k = n + 1, respectively. Consequently, 

{— sinh(fci^), at 1 < k < n 

A |sinh[(iV + 1 - k)tp] - coth (j^J cosh[(iV + 1 - k)ip] } , at n+l<k<N 

The mode factor ip and the amplitude A are determined by the conditions eq. (|3.32[) and eq. (|3.33p . The substitution 
of eq. (jXMl) in eqs. (|3~52l and P~53")) yields 

sinh[(n — l)ip] + aA |sinh[(iV — n)ip] — coth C^j cosh[(iV — n)</j] j + [1 + a — 2cosh(ip)] smh(nip) = 

asmh(nip) + A |sinh[(iV — n — l)tp] — coth (^J cosh[(iV — n — l)y>]| 

+ A[l + a - 2cosh(<^)] |sinh[(iV - n)ip] - coth cosh[(7V - n)<p]j = (3.40) 

At n = the first equation in eq. (|3.40[) becomes redundant and the second one can be written as 

|sinh[(A - l)ip] - coth cosh[(A - 1)^]| + [1 + a - 2cosh</?] jsinh(i\V) - coth (|) cosh(AV)} = (3.41) 

For the values 1 < n < N — 1 the amplitude A may be excluded from eq. (|3.40|) which leads to the equation 

a 2 smh(nip) = {sinh[(n — l)(p] + [1 + a — 2cosh</9] sinh(n<p)} 

f smhffiV - n - l)tp] - coth(w/2) coshffiV - n - l)ip] , , .] 

x <^ ^ — - ^ % _iri + i + a _2coshwH 3.42 

\ sinh[(AT - n)tp] - coth(<^/2) cosh[(7V - n)(p] 1 ^ J J V ' 

The substitution of n = in eq. f|3.42[) gives back eq. Q3.4ip as required by consistency. 

From the solution of the transcendental eq. (|3.41l) (for n = 0) or eq. (!3.42p for 1 < n < N — 1 one obtains the mode 
factor ip as a function of n and N, i.e. ip{N,n). Knowing ip, one may calculate the negative eigenvalue given by 
eq. (l3.35p . Making use of this in the equation for the factor k (see eq. (|2.28[) ) one obtains 



K(N,n) = l 



/ 80 2 
l-Wl + -^[cosh ¥ j(JV,n)-l] 



< (3.43) 



Finally, the first equation in (I3.40[) makes it possible to calculate the amplitude A as a function of N and n, i.e. 

A(N n) - sinh [(" ~ + I 1 + a Z 2 cosh (y)] smh{rvp) 
[ ' ' a{sinh[(A-n)^] -coth(^/2)cosh[(A-n)(^]} k ' ' 

As mentioned in Sec. MI CI the chain breaks when at least one bond approaches the position of the potential 
maximum, i.e. becomes "endangered" . On the other hand, within the Kramers approach this is a relatively rare 
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event, so that a simultaneous occurrence of a second, third, etc. "endangered" bonds may be neglected. Therefore, 
by calculating the rate constant for the total chain one should average over all possible locations of an "endangered" 
bond. Consequently, the rate constant of the total chain (see eq. (|2.21j) is given as 

1 1 N ^ 

k = _ e -PE bR _Y j \ K{N ,n)\ (3.45) 

where R and K(N,n) are given by eqs. p.30[) and p.43|) . respectively. By means of a MD-simulation, one may 
calculate the average first passage time r for crossing the barrier which is related in turn as k = (2r) _1 to the rate 
constant (see Sec. Ill CI and ref. |35j|V 

IV. SIMULATION RESULTS 

We present here the results our extensive MD simulations in order to verify the theoretical predictions and explore 
the limitations of the analytical treatment. Energy is given in units of D and length is measured in units of 1/a, 
the parameters of the Morse potential U(y) = D(l — e~ ay ) 2 . Mass is measured in terms of m, the mass of the bead. 
Time is measured in units of a^ 1 y / rn/D; temperature is measured in units of D/ks where the Boltzmann constant, 
ks has been set equal to 1.0. Unless otherwise mentioned, the ratio of the barrier height to temperature Eb/T has 
been set to 5, the value of the externally applied pulling force has been set equal to 0.25 and the friction coefficient 
of the Langevin thermostat, used for equilibration, has been set at 0.25. The integration step is 0.002. 




FIG. 4: Snapshots of a chain with 30 beads fixed at the left end while the right end is pulled by a constant force: (a) an 
equilibrated initial conformation, (b) a broken chain with the beads at the scission site shown in white. 

We start the simulation with each bead separated by a distance equal to the equilibrium separation of the effective 
bond potential U{y n ) (eq. (|3.2p h we then let the chain equilibrate with its environment using a Langevin thermostat. 
The number of integration steps for equilibration of a chain with 10 atoms in 1-D is 20000 and the number of 
equilibrating integration steps is increased linearly as the chain length is increased. Because of the presence of the 
external pulling force, the stretched chain attains equilibrium only locally and not globally, i.e., it never turns into a 
coil. This justifies the linear increase in the number of equilibration steps with increase in chain length as opposed to 
"quadratic" as in the Rouse model. Once equilibration is achieved, time is set to zero and one measures the elapsed 
time (in MD time units) before any of the bond lengths extends beyond the distance separating the metastablc 
minimum from the maximum, i.e., until one of the beads crosses the barrier (see Fig. [4]). We repeat the above 
procedure for a large number of events (5 x 10 4 -j- 2 x 10 7 ) so as to sample the stochastic nature of rupture and 
calculate properties like the mean rate of rupture, the distribution of breaking bonds regarding their position in 
the chain, the (First Passage Time) FPT distribution, etc., for chains of different length in both 3d and Id. As 
one of our principle objectives, we also investigate the issue of chain recombination (self-healing) which is frequently 
observed after a scission event occurs. We demonstrate that the mere barrier crossing is not a reliable criterion for 
chain breakage since the majority of broken bonds are observed to recombine. Therefore we develop an unambiguous 
criterion for true rupture as illustrated in a later subsection. 

A. Chain Scission - Simulation Results 

We compare here the simulation results with our theoretical prediction for the rupture probability of the n-th bond 
in a chain with N bonds. The probability for bond scission (or, exit probability, in the more general context of Section 
HJ), is given by eq. ([2301) . 

In Fig. [5] the normalized rupture probability for chains (with TV — 10 and N = 30) is shown with respect to the 
consecutive number of the individual bonds. The theoretical prediction, which follows from the numerical solution of 
eqs. (|3.41[) . (I3.42|) and (|3. 431) . is given in the inset. Both the theory- and MD-results indicate that the pulled end of 
the chain and the bonds in its vicinity break more frequently due to more freedom than those around the fixed end. 
Generally, the probability of rupture decreases steadily from the pulled end to the fixed end. For the longer chain, 
the end effects are not felt by the middle part of the chain and the probability of rupture P(N, n) is nearly uniform 
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Normalized Bond Number 



Bond Number 



FIG. 5: Normalized rupture probability vs consecutive bond number for chains with length N, subject to tensile force / = 0.25, 
and friction 7 = 0.25: (a) ID results for Eb/ksT = 5 and N — 10, 30. The consecutive number of the bonds is normalized 
as n/N for convenience. Insert shows the theoretical prediction, (b) Theoretical prediction for N — 30 along with simulation 
results for the Morse potential and the Double Harmonic potential. Inset shows the same for a chain with Morse interactions 
in both ID and 3D. 



forming a plateau-like region all over the length of the chain except at the ends. This feature is more pronounced in 
the theoretical rather than in the MD results. Such a comparison between the theory and MD-simulation findings has 
not been done before and could be viewed as a test for the merits and shortcomings of Kramers approach. Moreover, 
we believe that this reflects the incorporation of the collective unstable mode given in Sec. MI Dl 

One may assume that the detected discrepancy between theory and MD-simulation can be ascribed to the use of 
harmonic approximation around the metastable minimum and unstable maximum of the effective bond potential. 
To show this, we replaced the effective potential by a Double Harmonic potential - a parabola, and an inverted 
parabola, which approximates the effective potential in shape. The comparison of the three results - theoretical 
prediction, simulation with Morse potential, and simulation with the Double Harmonic potential is shown in Fig. [SJd. 
Evidently, the rupture probability distribution for the Double Harmonic potential exposes also a plateau-like region 
and matches the theoretical prediction quite closely. So we infer that the absence of the plateau in case of the Morse 
potential can be ascribed to its anharmonicity. In the inset of Fig. [SJd, we show the simulation results with the Morse 
potential for a 30-atom chain in ID and 3D. Interestingly, one may see that the MD-results in 3D indicate a flatter 
distribution P(N,n)) than in ID, coming thus closer to theoretical predictions. As far as in 3D there exist much 
more configurations with endangered bonds due to transversal displacements of the beads, bond length fluctuations 
are suppressed as compared to 2D and the bond anharmonicity is less pronounced. 

Theory predicts that the first passage time distribution goes asymptotically as W ~ exp(— t/r) ( see eq. (|2.33|) ). In 
Fig. [SI we plot the FPT distribution for a chain with 30 beads. One should note here the considerable difference of 
W(t) between ID and 3D. The long time tail of the distribution is indeed seen to decay exponentially. We estimate 
the mean FPT (for the whole chain and not of an individual bond) as t\d = 83.2 and t^d = 36.0. The theoretical 
estimate for a 30-bead chain in ID is t\d — 866.5 which is an order-of-magnitude larger than that estimated from 
the simulations. This finding is in aj e I in (26[. 




FIG. 6: First passage time distributions W(t) against elapsed time r in ID and 3D for a chain with N = 30, / = 0.25 and 
7 = 0.25. The inset shows the same in semilog coordinates. 
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While assessing these results one should bear in mind that, as mentioned above, an event of bond scission can 
be defined as " barrier crossing" whereby an endangered bond stretches beyond the position of the maximum of our 
potential U(y). As has been discussed in the previous section, this position of the barrier depends on the externally 
applied force once the parameters of the Morse potential are fixed. For an applied force of magnitude 0.25, the barrier 
position corresponds to a critical bead-bead separation (bond length) of 2.92. Since frequently such scission event is 
immediately followed by recombination (that is, by self-healing of the bond), in what follows we have also considered 
a more stringent criterion for irreversible rupture. As a possible choice for the critical bead-bead separation one may 
take the larger value of 5 which renders self-healing events virtually improbable (see below) . 

It is of interest to examine the dependence of mean FPT (r) on temperature T for a given applied force /. In FigJT] 
we show the change in FPT for irreversible rupture (and not barrier crossing) on inverse temperature for / = 0.25 
and 7 = 0.25. As expected from eq. (|3.45p . the Arrhenian nature of this relationship is clearly manifested in a semilog 
plot. From the exponentially fitted values, the effective barrier height is estimated to be 0.176 in ID and 0.205 in 3D. 
The theoretical estimate of the barrier for a force of magnitude 0.25 is AEb ~ 0.27. The somewhat lower value of the 
effective barrier has also been observed earlier in simulations with the Lennard-Jones potential for the fixed strain 
ensemble (2lj . 
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FIG. 7: Mean first passage time (r) vs inverse temperature for a chain with N = 30, / = 0.25 and 7 = 0.25. 



Eventually, in Fig. [5] we show the iV-dependence of the mean FPT (r) and the average scission rate for both the 
Morse chain and the double-harmonic chain. We recall that we considered the unstable collective mode in Section 
MID as the one responsible for the rupture event. Therefore, one might expect that the TV-dependence of the scission 
rate would clearly reveal the degree of collectivity in bond breaking. Indeed, a linear increase in the scission rate 
with growing chain length would indicate the independent character of scission events, that is, the breakage of a single 
bond occurs irrespective of the neighboring bonds. In contrast, if the nature of bond scission is a strongly collective 
process, this iV-dependence should be negligible. Evidently, Fig. [8] is a clear manifestation of the latter, contrary to 
an earlier assumption [3, [IH, HH that the total probability for scission of polymer with N bonds is N times that 
of a single bond. One sees a negligible decline in the mean FPT in Fig. [5^, which appears as a very weak rise in 
the total rate in Fig. [8]d, much weaker (with slope ~ 0.1) than the presumed linear growth with N. A visible finite 
size effect is detected only for N — 10, indicating a strong influence of the boundaries (i.e., of the chain ends). This 
effect is easily understood from the inset to Fig. [5ji where the absence of a plateau suggests that the bonds in the 
short chain are not equivalent. One can also verify from Fig. [5] that the double-harmonic chain breaks faster than the 
Morse one, 



B. The breaking of a bond 



We now examine the expansion rate of the breaking bond. Theory predicts that the length of the endangered bond 
extends with time as exp(— ret) (see eq. (|2.14|l ). where re, the transmission factor given by eq. (|3.43[l . is negative. 

Exponential fits of the growth curves around the position of the barriers (~ 2.92 for / = 0.25 and ~ 3.05 for 
/ = 0.15) give |re| ~ 0.01 for / = 0.15 and |re| ~ 0.02 for / = 0.25 is given in FigJH Analytical calculations yield 0.11 
and 0.22 - Fig. |9]- respectively, an order of magnitude higher than that estimated from the simulations and consistent 
with the result obtained earlier from the lifetime distribution. 
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FIG. 8: Mean first passage time (r) (a) and average scission rate (b) vs chain length N for Morse and double-harmonic chains 
with / = 0.25 and 7 = 0.25. 




FIG. 9: Stretching of a breaking (15-th) bond with time in a chain with 30-beads for / = 0.15 and / = 0.25: (a) Evolution 
of the lengths of a breaking bond and of its neighbors. The dashed line denotes an extension equal to the barrier position, 
(b) Variation of the length of a breaking bond with time in semi-logarithmic coordinates. Dashed lines indicate exponential 
increase in agreement with theoretical predictions. The inset shows the expansion in normal coordinates for visual aid. The 
length growth has been monitored up to an expansion of 15, well beyond the location of the barrier. Data has been averaged 
over many events. The dashed black line shows the position of the barrier while the solid line marks an expansion of 5.0, later 
set as the criterion for irreversible rupture for an applied force / = 0.25. 



Eventually, we show the rate of expansion of the monitored bond for two different forces, note that the exponential 
character of the curves is best pronounced when the growing bond length is in the vicinity of the barrier (hump) 
position which is indicated in Fig. |H]by horizontal lines, in agreement with the theoretical predictions, eq. (j2.14j) . The 
speed of expansion itself is displayed in Fig. [10] (statistic fluctuations are strongly pronounced for the weaker tension 
/ = 0.15) which demonstrates that also the speed of bond extension attains a maximum at the moment it goes over 
the barrier position. 
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FIG. 10: Mean expansion rate of the breaking 15-th bond. 



V. SELF-HEALING (RECOMBINATION) OF BROKEN BONDS 



The process of bond breakage is not fully described by the escape from the metastable well and the barrier crossing. 
Once the barrier has been surmounted, the " endangered bond" could further increase its length up to a critical value 
which still leaves the possibility for eventual return and healing of the chain. In contrast, beyond this critical value 
the healing probability gets negligibly small and the chain breaks irreversibly. This healing process has been pointed 
out earlier [20j and observed in a MD-simulation Below we suggest a theoretical description of the process of 
recombination and derive analytical expressions for the healing time and distance distribution function. One should 
bear in mind that it is the thermal fluctuations that initiate recombination events while the downhill ramp potential of 
the applied tensile force always acts in the opposite "destructive" direction. This potential may still be approximated 
by an inverted parabola before the " endangered bond" length approaches a critical size whereby the two pieces of the 
chain would move apart irreversibly. 



A. Theory 



As far as the healing process is largely localized on the " endangered bond" , it may be treated as an effectively one- 
dimensional problem. The corresponding Fokker-Planck equation for the probability distribution function of positions 
(in our case - "endangered" bond lengths) and velocities (bond length velocities) is known as Kramers equation f2o]] 
and has been used to describe reaction kinetics. The complete solution of the Kramers equation for the harmonic as 



well as inverted harmonic (or inverted parabolic) potential was Fig. [8] given by Risken 31 [ 

Consider the Kramers equation for the inverted harmonic potential U(x) — — Sl|a; 2 /2, assumed to represent the top 
of the barrier. The coordinate-velocity transition probability P(x,v,t\x' ,v' ,0) is governed by the following equation 
(see Sec. 10 in (U): 



2 



§ i P = -l [v P ] + l m >( X ) + , V )P ]+ , V l^P (5.1) 



where the thermal velocity i>th = y/T/m and the initial conditions are fixed as P{x, v, t — Ola;', v' , 0) = S(x—x')S(v—v'). 
The general solution of eq. (|5.1|) is given by [3l| and is relegated to Appendix O 



1. Healing time distribution function 



One may treat the healing time distribution function by noting that the process starts at the top of the inverted 
harmonic potential U(x) = —VB^x 1 j1^ where x' — and v' — \k\ . The coordinate of the turning point (healing 
distance) is counted with respect to x — whereas the velocity at that point might be arbitrary within the interval 
— oo < v < (i.e., the return velocity is pointing to the left). Thus, the healing time PDF can be obtained from the 
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transition probability eq. (IC1[) as 

o 

fUi(*)= / d«P(0,«,t|0,|/s|,0) 

— oo 

After explicit integration in eq. (|5.2[) and taking into account eq. (|C6|) . one obtains 



(5.2) 



1 



iW(t) = o 



1 



2 W 27r(T ;ca ;(t) 



exp 



2a xx (t) 



1-erf 



(t) v(t) - <J xv (t) x(t) 

\]2o xx (t) deter 



(5.3) 



where the error function erf(z) = (2/^/n) J* e x dx, and £c(i), v(i) are given by eq. (|C7|) (with x' = and w' = |/c| 



2. Healing distance distribution function 

In order to determine the distribution of distances h from the hump at x — where the extending bond may still 
turn back to shrinking, one may consider the PDF for the maximal divergence distance Q(h) before the bond returns 
and heals. In this case a critical value h c can be defined as a point where Q(h) gets very small. 

One may express Q(h) in terms of the transition probability P(x,v,t\x' ,v' ,t') by taking into consideration the 
following arguments. The overall return process can be seen as the composition (recall that the process is Markovian) 
of two processes. The first one starts from the top of the inverted harmonic potential with the initial x 1 = and 
v' = \k\ and continues down the ramp until a turning point x = h where the velocity becomes zero, i.e., v = 0, at an 
intermediate time moment t 1 . The probability of this process is given by P(h, 0, t'\0, \k\,0). The reverse process starts 
from the turning point (i.e., x 1 = h,v' = 0) at time moment t' , and continues back to the top of the potential where 
x = and the velocity can be any in the interval — oo < v < 0. The latter means that the corresponding transition 
probability should be integrated over velocity, i.e., f_ dvP(0, v, t\h, 0, t'). Finally, since the total time interval may 
also be arbitrary, one should integrate over times. Thus, one obtains 



oo t u 

Q(h) = J dt J dt' J dv P(0,v,t\h,0,t') P(h,0,t'\0, \k\,0) 



-oo 
oo 



= / dr I dvP(0,v,T\h,0,0) / dt'P(h,0,t'\0, \k\,0) = 1(h) JQi) 



i(h) 



j(h) 



where we have changed the order of time integration. Taking into account the form of the transition probability, eq. 
(|C1[) . and after integrating over the velocity, the expression for 1(h) reads 



1(h) = I I dtj - 1 . . 
v ; 2 J V 2 ™ xx (t) 



exp 



2<J xx (t) 



1 — erf 



(t)v(t) - a xv (t)x(t) 
\j2a xx (t) deter 



(5.4) 



where x(t) and v(t) are given by 

x(t) = G xx (t)h 

v(t) = G vx (t)h (5.5) 
By making use of the expressions for the inverse cr-matrix, eq. (|C6[) . the expression for J(h) takes on the form 



J(h) 



1 

2^ 



dt- 



1 



v / det( 



exp 



with 



a vv (t)[h - x(t)] 2 + 2a xv (t)[h - x(t)]v(t) + a xx (t)v 2 (t) 
2 deter 



x(t) = G xv (t)\n\ 
v(t) = G vv (t)\n\ 



(5.6) 



(5.7) 
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B. MD-simulation of Self-Healing 

The simulation scheme for sampling the self-healing distributions is as follows - once the largest bond in the chain 
crosses the barrier position, one monitors its expansion for 1000 integration steps and records the maximum expansion 
and the time spent before it crosses back the barrier. Frequently such a pseudo-broken bond heals again so that from 
the above record, one may recover the distribution of the maximum expansion beyond the barrier, and the time spent 
in a pseudo-broken state before the healing occurs. 




12 3 4 



Expansion 

FIG. 11: Probability distribution of healing bond lengths and healing times (inset) for a 30-bead chain in ID and / = 0.25. 
Both plots are shown in semilog coordinates. Dashed lines denote theoretical predictions. 

In Fig. [TTJwe show the distribution of the healing lengths and healing times (inset) obtained from simulations and 
from theory. Evidently, the healing time distribution Pkeax{t) looks qualitatively identical for both the theoretical 
treatment (which is based on eq. (|5.3[) )and the simulations: there is a maximum just beyond the barrier crossing 
event for t =/= and a fast (exponential) decay thereafter. The existence of the maximum can be attributed to 
the thermostat - just after crossing the barrier one has to wait for a while a thermal kick turns the trajectory into 
reverse direction. On the other hand, for longer time intervals the downhill motion wins (i.e., the healing becomes 
progressively improbable) and Pheai(^) decreases rapidly with time. 
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FIG. 12: The healing expansion distribution Q(h) for / = 0.15 and / = 0.25 plotted in semi-logarithmic coordinates. The zero 
of the x-axis corresponds to the barrier positions, which are different for the two different values of the force. 

The distribution of the healable expansions beyond the barrier is of greater importance because it sets a more 
stringent criterion of true (irreversible) rupture. We see that for / = 0.25, the probability of healing is negligible for 
an expansion of ~ 2.0 beyond the barrier, i.e., a bond expansion of 2.92 + 2.0 « 5.0. Once the expansion of the bond 
reaches h c ~ 5.0 for / = 0.25, it is highly unlikely for the bond to heal. One should note that the process is stochastic 
in nature so that a bond that has expanded even beyond h c may still manage to heal. Such an event, however, is 
highly unlikely and the criterion defined in the above manner serves to be a good one for all practical purposes. The 
theoretical distribution (calculated according to eqs. (|5.4p . (|5.4p and (|5.6|1 ). decays much faster than that obtained 
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from simulations; this is because of the faster descent of the quadratic potential (considered in theory within the 
linearized approximation) as compared to the roughly linearly decaying effective potential beyond the barrier. 

How does the above picture change when the conditions of the simulation are varied? To answer this question, we 
have performed simulations and we find that the above distribution does not vary much with friction, chain length or 
whether the simulation is performed in ID or 3D. The distribution however depends on the applied tensile force as 
shown in Fig. [12] for two different forces. 

As expected intuitively, the distribution decays faster for a larger force because healing is more unlikely in case of a 
steeper unstable potential beyond the barrier. From a practical point of view, such a distribution for several different 
forces can be very helpful in determining a criterion for irreversible rupture. 

We next turn to the question of the dependence of the healing on temperature. 



N=10 f=0.25 7=0.25 




FIG. 13: (a) The healing fraction for / = 0.25 and 7 = 0.25 at different values of the temperature. A power-law regression 
gives an exponent of ~ 0.2 in ID and ~ 0.3 in 3D. (b) Increase of the fraction of self-healing bonds with friction 7. 

Figure [T31 shows the fraction of self-healing events on the temperature for / = 0.25 and 7 = 0.25. As expected, 
healing becomes more frequent as the temperature is increased because the bead that has crossed the barrier now 
receives stronger thermal kicks capable of sending it back over the barrier. One can also see that for the same 
temperature, healing is (understandably) more frequent in II? than in 31?. The effect of growing friction 7 is similar, 
on the one hand it strengthens kicks by the thermostat on the beads at the verge of breaking, and on the other hand 
it delays such beads in "rolling down" the downhill ramp potential, providing more time to receive a thermal kick 
which could turn them back. 



VI. CONCLUSION 

In the present work we have studied the process of polymer chain breakage for linear chains subjected to constant 
tensile force by multidimensional Kramers-Langer approach and compared theoretical predictions to results of exten- 
sive MD simulations in one- and three dimensions. The adopted theoretical treatment makes it possible to consider 
collective unstable modes as being mainly responsible for chain scission. Comparison with simulation data as, for 
example, the distribution of the probability for scission over bond index, the scission time distribution, variation of 
MFPT with temperature and chain length as well as the bond breaking dynamics indicates that the Kramers-Langer 
approach agrees qualitatively with observations. 

We demonstrate that the recombination (self-healing) dynamics of the breaking bonds can be qualitatively repro- 
duced by the Kramers equation. We derive analytic expressions for the healing time and expansion distributions in 
reasonable agreement with the MD data and reveal the variation of the fraction of healing bonds with temperature 
and friction. In fact, we demonstrate that more that 50% of the crossings of the activation barrier end up as healed 
bonds again (this fraction changes with temperature) especially in the 3D simulations. This finding should be kept 
in mind in the assessment of earlier research work where only ID simulations of breaking chains have been performed 
and bond-healing was not allowed for. 

One should emphasize, however, that we still find a discrepancy of nearly an order of magnitude between theory 
and computer experiment regarding the mean rate of scission, that is, the rate of bond breaking is significantly 
underestimated by theory (cf . also [26| ) . We believe that the origin of this discrepancy is related to the existence of non- 
linear localized excitations (breathers) in the anharmonic lattice which remain out of the scope of this investigation. 
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Recent studies [391 ] indicate the breathers play very important role regarding energy transfer in discrete non-linear 
chains which serve a generic model for biopolymers. Clearly, furter research in this direction is needed before a good 
understanding of the whole problem is reached. 
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Appendix A: Calculation of ^\ ViTijVj 



Let us calculate a = J2i j viTijVj which occurs in eq. (I2.44[) . Because J2i=i v 1 = 1 this can be seen as a eigenvalue 
problem, i.e. 



2iV 



(Al) 



The corresponding characteristic equation reads 

d 

The matrix Tij has a block structure (see eq. (|2.23[1 ) so that eq. (| A2|) takes the form 



dct[Fy — 5ij<j] — 
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(A3) 



On the other hand the determinant of he block-diagonal matrix det ( ^ ^ ) = det A det D = (—<r) N (mj — a) N = 0. 
As a result 



a = ffl/y 



(A4) 



Appendix B: How to calculate the determinant of a symmetrical tridiagonal matrix 

The typical tridiagonal N x iV-matrix which is common in the context of one-dimensional string of beads model 
has the following form 



T(N) = 



( 2 1 

-1 2 

-1 



V 



0-12 -1 
-1 !-/?/ 



(Bl) 



where f3 is an arbitrary rational number. To calculate the determinant we use the Laplace's formula (3 31 ] which holds 
that for an arbitrary matrix A the determinant det(A) = Ylj=i Aij(—iy +J Mij where is the minor (i.e. the 
determinant of the matrix that results from A by removing the i-th row and the j-th column) . By making use this 
formula for the matrix in eq. (IB1|) with respect to the first row and expanding similarly again, one finds the recurrence 
relation 



det[T(iV)] = 2dct[T(A r - 1)] - dct[T(/V - 2)] 



(B2) 
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This reccurence relation to be solved needs two initial conditions 

det[T(l)] = l~/3 det[T(2)] = det ^ _J ^"^=1-2/3 (B3) 

Thus 

det[T(3)] = 2(1-2/3) - 1 + /3 =1-3/3 
det[T(4)] = 2(1-3/3) - 1 + 2/3 =1-4/3 

(B4) 



det[T(AT)] = 1-/3 AT 
As a result we can finally write down 

det[T(A0] = (1 - 0N) (B5) 

Appendix C: Solution of one-dimensional Kramers equation 

The general solution of eq. (|5.1[> reads [3l[ 

P(a^,t|z'y,0) = ^^172 ^{-jl^Wla^-^f"^^..!*-^)] [«-«(*)] 

- ^-^Wb-^)] 2 } (ci) 

where the 2x2 c-matrix 



has the following elements has the following elements 

<?xx(t) 



"< v th [A1 + A2 i 4 c-( Al+A2 )* - 1 e~ 2Al * e~ 2Aat 



(Ai - X2) 1 



A1A2 Ai + A2 Ai 



*™(t) = tt 7% ; , 2 (Ai + A 2 + i^f- [ e -^+ A = - 1] - Ai e- 2A ^ - A 2 e" 2A 4 (C3) 
(Ai-A 2 ) 2 ( A1+A2 J 



and where the eigenvalues Ai and A 2 are 

1 



Ai = ^( 7 +x/7 2 +4tt 2 ) (C4) 



A 2 = -( 7 -^7 2 +4^) = ^<0 (C5) 

We underline that the negative eigenvalue Ai is nothing but the transmission factor k given by eq. ()3.43|) whereas 
the characteristic frequency Q 2 is given by eq. (|3.5[) . The determinant and the inverse a- matrix in eq. (|C1|) arc 
defined as as follows 

deter = o xx <Jw — cr xv 
(cr 1 )xx 



det rr 

(f )xv = (c 1 )vx = 



det rr 



= (C6) 
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The mean values x(t) and v{t) in eq. (IC1[) are given by 



x(t) 
v(t) 



G xx {t) x' + G xv (t) v' 
G vx (t) x' + G vv (t) «' 



(C7) 



where the Green function matrix elements read 



G X x(t) 



G X v(t) 
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